rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)


#This is to obtain:
#-figure_1a
#-figure_1b
#-figure_1c
#-figure_1d
#-figure_a6a
#-figure_a6b
#-figure_a7a
#-figure_a7b
#-figure_a7c
#-figure_a7d
#-figure_a8a
#-figure_a8b
#-figure_a8c
#-figure_a8d
#-figure_a9a


library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
#library("rgdal")
library("corrplot")
library("Hmisc")
library("glue")
library("ggplot2")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
#library("rgeos")
library("ggrepel")
library("lemon")
library("sandwich")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("pgirmess")
library(sf)
library("tidyverse")
library("geojsonsf")


data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")


#Step1. Turning some variables to numeric
data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)


data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_total<-data_cov_mob$gspo+ data_cov_mob$gkul+ data_cov_mob$gnat
data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
summary(data_cov_mob$bridging_tot_pop)
mean_brid<-summary(data_cov_mob$bridging_tot_pop)[4]
data_cov_mob$bridging_tot_pop_dummy<-ifelse(data_cov_mob$bridging_tot_pop>mean_brid, 1, 0)

#Bonding Networks
data_cov_mob$bonding_total<-data_cov_mob$gsoz+ data_cov_mob$gpol+ data_cov_mob$ginter
data_cov_mob$bonding_tot_pop<-data_cov_mob$bonding_total/data_cov_mob$pop_total*1000
summary(data_cov_mob$bonding_tot_pop)
mean_bond<-summary(data_cov_mob$bonding_tot_pop)[4]
data_cov_mob$bonding_tot_pop_dummy<-ifelse(data_cov_mob$bonding_tot_pop>mean_bond, 1, 0)

#All Networks
data_cov_mob$gvereine<-as.numeric(data_cov_mob$gvereine)
data_cov_mob$vereine_tot_pop<-data_cov_mob$gvereine/data_cov_mob$pop_total*1000


summary(data_cov_mob$vereine_tot_pop)
mean_ver<-summary(data_cov_mob$vereine_tot_pop)[4]
data_cov_mob$vereine_tot_pop_dummy<-ifelse(data_cov_mob$vereine_tot_pop>mean_ver, 1, 0)
table(data_cov_mob$vereine_tot_pop_dummy)
table(data_cov_mob$vereine_tot_pop)

data_cov_mob_week<-subset(data_cov_mob, week=="01")
data_cov_mob_week<-subset(data_cov_mob_week, select = c("GEN", 
                                                        "vereine_tot_pop",
                                                        "bridging_tot_pop", "bonding_tot_pop"))

data<-st_drop_geometry(data_cov_mob)
data$week_new<-as.numeric(data$week)+1
data$week_new<-paste0("0", data$week_new)  
data$week_new2<-ifelse(nchar(data$week_new)==3, str_remove(data$week_new, "^0+"), data$week_new)
data$week<-data$week_new2
data<-subset(data, select = -c(week_new, week_new2))
data$post_treat<-ifelse(as.numeric(data$week)>10, 1,0)


data$post_treat_vereine_tot_pop_dummy<-data$post_treat*data$vereine_tot_pop_dummy
data$post_treat_bridging_tot_pop_dummy<-data$post_treat*data$bridging_tot_pop_dummy
data$post_treat_bonding_tot_pop_dummy<-data$post_treat*data$bonding_tot_pop_dummy
data$post_treat_pct_service_sec_narrow<-data$post_treat*data$pct_service_sec_narrow

unique(data$post_treat_vereine_tot_pop_dummy)
unique(data$post_treat_bridging_tot_pop_dummy)
unique(data$post_treat_bonding_tot_pop_dummy)
unique(data$post_treat_pct_service_sec_narrow)


table(data$post_treat, data$mean_mobil)


df <-data%>%
  group_by(post_treat) %>%
  dplyr::summarise_at(vars(mean_mobil), list(name = mean))


#############################
#Function for calculating SE#
#############################

# Calculate robust confidence intervals
se_robust_cluster <- function(x){
  res<-coeftest(x, vcov = cluster.vcov(x, data$kreis_name, force_posdef=T))[, 2]
  return(res)}

ci_robust_cluster <- function(x){
  res_ci<-coefci(x, vcov = cluster.vcov(x, data$kreis_name, force_posdef=T))
  return(res_ci)}

the_FE <-names(data[ , grepl( "SN_L_" , names(data))])
the_FE<-the_FE[-length(the_FE)]
small_FE<-c('SN_L_01',
            'SN_L_02',
            'SN_L_03',
            'SN_L_05',
            'SN_L_06',
            'SN_L_12',
            'SN_L_13',
            'SN_L_14',
            'SN_L_15')


######################
#Bridging Time Trends#
######################

df2 <- subset(data, select = c(mean_mobil, week, 
                               bridging_tot_pop_dummy))
df2<-subset(df2, !is.na(bridging_tot_pop_dummy))
df2$type<-ifelse(df2$bridging_tot_pop_dummy==1, "Above Mean", "Below Mean")
unique(df2$type)

df2<-st_drop_geometry(df2)
df3 <- df2%>%
  group_by(type, week)%>%
  summarise(mean_mobil = mean(mean_mobil))


df2_covid <- subset(data, select = c(covid_diff, week, 
                               bridging_tot_pop_dummy))
df2_covid<-subset(df2_covid, !is.na(bridging_tot_pop_dummy))
df2_covid$type<-ifelse(df2_covid$bridging_tot_pop_dummy==1, "Above Mean", "Below Mean")
unique(df2_covid$type)

df2_covid<-st_drop_geometry(df2_covid)
df3_covid <- df2_covid%>%
  group_by(type, week)%>%
  summarise(covid_diff = mean(covid_diff))


x<-seq(-10, 42, by = 1)
length(x)
df3$observation<-NA
df3$observation[df3$type=="Above Mean"] <- x
df3$observation[df3$type=="Below Mean"] <- x

df3_covid$observation<-NA
df3_covid$observation[df3_covid$type=="Above Mean"] <- x
df3_covid$observation[df3_covid$type=="Below Mean"] <- x

df3$new_type<-NA
df3$new_type[as.numeric(df3$week)<11 & df3$type == "Above Mean"]<-"Above Mean Pre"
df3$new_type[as.numeric(df3$week)<11 & df3$type == "Below Mean"]<-"Below Mean Pre"
df3$new_type[as.numeric(df3$week)>=11 & df3$type == "Above Mean"]<-"Above Mean Post"
df3$new_type[as.numeric(df3$week)>=11 & df3$type == "Below Mean"]<-"Below Mean Post"


df3_covid$new_type<-NA
df3_covid$new_type[as.numeric(df3_covid$week)<11 & df3_covid$type == "Above Mean"]<-"Above Mean Pre"
df3_covid$new_type[as.numeric(df3_covid$week)<11 & df3_covid$type == "Below Mean"]<-"Below Mean Pre"
df3_covid$new_type[as.numeric(df3_covid$week)>=11 & df3_covid$type == "Above Mean"]<-"Above Mean Post"
df3_covid$new_type[as.numeric(df3_covid$week)>=11 & df3_covid$type == "Below Mean"]<-"Below Mean Post"



df3$week
as.numeric(df3$week)
df3_small<-subset(df3, as.numeric(df3$week)>=5 & as.numeric(df3$week)<=16)
irislabs1<-seq(min(as.numeric(df3_small$week)), max(as.numeric(df3_small$week),1))
#irislabs2<-unique(df3_small[as.numeric(df3_small$week) %in% irislabs1,]$t)


df3_small_covid<-subset(df3, as.numeric(df3_covid$week)>=5 & as.numeric(df3_covid$week)<=16)
irislabs1_covid<-seq(min(as.numeric(df3_small_covid$week)), max(as.numeric(df3_small_covid$week),1))

df3$week<-as.numeric(df3$week)
sample_abom_bel10<-subset(df3, week<=10 & type=="Above Mean")
sample_abom_abo10<-subset(df3, week>10 & type=="Above Mean")
sample_belm_bel10<-subset(df3, week<=10 & type=="Below Mean")
sample_belm_abo10<-subset(df3, week>10 & type=="Below Mean")


df3_covid$week<-as.numeric(df3_covid$week)
sample_abom_bel10_covid<-subset(df3_covid, week<=10 & type=="Above Mean")
sample_abom_abo10_covid<-subset(df3_covid, week>=10 & type=="Above Mean")
sample_belm_bel10_covid<-subset(df3_covid, week<=10 & type=="Below Mean")
sample_belm_abo10_covid<-subset(df3_covid, week>=10 & type=="Below Mean")


mo_abom_bel10 <- lm(mean_mobil ~ week, data=sample_abom_bel10)
summary(mo_abom_bel10)
xseq<-seq(from=min(sample_abom_bel10$week), to=max(sample_abom_bel10$week))
pred <- predict(mo_abom_bel10)
df_abom_bel10 <- data.frame(week = xseq, pred)
df_abom_bel10$type_specific<-"Above Mean Pre"
df_abom_bel10$type<-"Above Mean"

mo_abom_bel10_covid <- lm(covid_diff ~ week, data=sample_abom_bel10_covid)
xseq_covid<-seq(from=min(sample_abom_bel10_covid$week), to=max(sample_abom_bel10_covid$week))
pred_covid <- predict(mo_abom_bel10_covid)
df_abom_bel10_covid <- data.frame(week = xseq_covid, pred_covid)
df_abom_bel10_covid$type_specific<-"Above Mean Pre"
df_abom_bel10_covid$type<-"Above Mean"

mo_abom_abo10 <- lm(mean_mobil ~ week, data=sample_abom_abo10)
xseq<-seq(from=min(sample_abom_abo10$week), to=max(sample_abom_abo10$week))
pred <- predict(mo_abom_abo10)
df_abom_abo10 <- data.frame(week = xseq, pred)
df_abom_abo10$type_specific<-"Above Mean Post"
df_abom_abo10$type<-"Above Mean"


mo_abom_abo10_covid <- lm(covid_diff ~ week, data=sample_abom_abo10_covid)
xseq_covid<-seq(from=min(sample_abom_abo10_covid$week), to=max(sample_abom_abo10_covid$week))
pred_covid <- predict(mo_abom_abo10_covid)
df_abom_abo10_covid <- data.frame(week = xseq_covid, pred_covid)
df_abom_abo10_covid$type_specific<-"Above Mean Post"
df_abom_abo10_covid$type<-"Above Mean"


mo_belm_bel10 <- lm(mean_mobil ~ week, data=sample_belm_bel10)
xseq<-seq(from=min(sample_belm_bel10$week), to=max(sample_belm_bel10$week))
pred <- predict(mo_belm_bel10)
df_belm_bel10 <- data.frame(week = xseq, pred)
df_belm_bel10$type_specific<-"Below Mean Pre"
df_belm_bel10$type<-"Below Mean"


mo_belm_bel10_covid <- lm(covid_diff ~ week, data=sample_belm_bel10_covid)
xseq_covid<-seq(from=min(sample_belm_bel10_covid$week), to=max(sample_belm_bel10_covid$week))
pred_covid <- predict(mo_belm_bel10_covid)
df_belm_bel10_covid <- data.frame(week = xseq_covid, pred_covid)
df_belm_bel10_covid$type_specific<-"Below Mean Pre"
df_belm_bel10_covid$type<-"Below Mean"


mo_belm_abo10 <- lm(mean_mobil ~ week, data=sample_belm_abo10)
xseq<-seq(from=min(sample_belm_abo10$week), to=max(sample_belm_abo10$week))
pred <- predict(mo_belm_abo10)
df_belm_abo10 <- data.frame(week = xseq, pred)
df_belm_abo10$type_specific<-"Below Mean Post"
df_belm_abo10$type<-"Below Mean"


mo_belm_abo10_covid <- lm(covid_diff ~ week, data=sample_belm_abo10_covid)
xseq_covid<-seq(from=min(sample_belm_abo10_covid$week), to=max(sample_belm_abo10_covid$week))
pred_covid <- predict(mo_belm_abo10_covid)
df_belm_abo10_covid <- data.frame(week = xseq_covid, pred_covid)
df_belm_abo10_covid$type_specific<-"Below Mean Post"
df_belm_abo10_covid$type<-"Below Mean"


new<-rbind(df_abom_bel10,
           df_abom_abo10,
           df_belm_bel10,
           df_belm_abo10)
new$line<-"predicted"
df3$line<-"actual"
names(new)[names(new)=="pred"]<-"mean_mobil"
names(new)[names(new)=="type_specific"]<-"new_type"

df3 <- df3[, c("week", "mean_mobil", "new_type", "type", "line")]
head(df3, 3)
head(new, 3)
df4<-rbind(df3, new)
x<-seq(-10, 42, by = 1)
length(x)
df4$observation<-NA
df4$observation[df4$type=="Above Mean"] <- x
df4$observation[df4$type=="Below Mean"] <- x


new_covid<-rbind(df_abom_bel10_covid,
                 df_abom_abo10_covid,
                 df_belm_bel10_covid,
                 df_belm_abo10_covid)
new_covid$line<-"predicted"
df3_covid$line<-"actual"
names(new_covid)[names(new_covid)=="pred_covid"]<-"covid_diff"
names(new_covid)[names(new_covid)=="type_specific"]<-"new_type"


df3_covid <- df3_covid[, c("week", "covid_diff", "new_type", "type", "line")]
head(df3_covid, 3)
head(new_covid, 3)
df4_covid<-rbind(df3_covid, new_covid)
x<-seq(-10, 42, by = 1)
length(x)
df4_covid$observation<-NA
df4_covid$observation[df4_covid$type=="Above Mean"] <- x
df4_covid$observation[df4_covid$type=="Below Mean"] <- x


df4$observation<-as.character(df4$observation)
df4_covid$observation<-as.character(df4_covid$observation)

selection<-df4$observation[!stringr::str_detect(df4$observation, '-')]
selection_covid<-df4_covid$observation[!stringr::str_detect(df4_covid$observation, '-')]


df4$observation[!stringr::str_detect(df4$observation, '-')]<-paste("+", selection, sep = "")
df4$t<-paste("t", df4$observation, sep = "")
df4$t[df4$t=="t+0"]<-"t"

df4_covid$observation[!stringr::str_detect(df4_covid$observation, '-')]<-paste("+", selection_covid, sep = "")
df4_covid$t<-paste("t", df4_covid$observation, sep = "")
df4_covid$t[df4_covid$t=="t+0"]<-"t"

####################
#Trend All Bridging#
####################
irislabs1<-seq(1,53,6)
irislabs2<-unique(df4[as.numeric(df4$week) %in% irislabs1,]$t)
length(irislabs2)
length(irislabs1)


trend_all_bri<-ggplot(df4, aes(x=as.numeric(week), y=mean_mobil, linetype = line, color=type))+
  geom_line()+
  scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
trend_all_bri<-reposition_legend(trend_all_bri, 'bottom left')
ggsave(trend_all_bri, file = "./paper/graphs/figure_1c.jpg", 
       height = 12, width = 20, 
       units = "cm", dpi = 200)



df4_covid_nonpredicted<-subset(df4_covid, line=="actual")
trend_all_bri_covid<-ggplot(df4_covid_nonpredicted, aes(x=as.numeric(week), y=covid_diff, linetype = line, color=type))+
  geom_line()+
  scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))

trend_all_bri_covid<-reposition_legend(trend_all_bri_covid, 'top left')
ggsave(trend_all_bri_covid, file = "./paper/graphs/figure_a6a.jpg", 
       height = 12, width = 20, 
       units = "cm", dpi = 200)


####################
#Trend Zoom Bridging#
####################
df3_small<-subset(df4, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df3_small[as.numeric(df3_small$week) %in% irislabs1,]$t)
length(irislabs2)

trend_zoom_bri<-ggplot(df3_small, aes(x=as.numeric(week), y=mean_mobil, linetype = line, color=type))+
  geom_line()+
  scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  scale_x_continuous(name = "Weeks", breaks=seq(5,16,1))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
trend_zoom_bri<-reposition_legend(trend_zoom_bri, 'bottom left')
ggsave(trend_zoom_bri, file = "./paper/graphs/figure_1a.jpg", 
       height = 12, width = 20, 
       units = "cm", dpi = 200)


#####################
#Bonding Time Trends#
#####################


df2 <- subset(data, select = c(mean_mobil, week, 
                               bonding_tot_pop_dummy))
df2<-subset(df2, !is.na(bonding_tot_pop_dummy))
df2$type<-ifelse(df2$bonding_tot_pop_dummy==1, "Above Mean", "Below Mean")
unique(df2$type)

df2<-st_drop_geometry(df2)
df3 <- df2%>%
  group_by(type, week)%>%
  summarise(mean_mobil = mean(mean_mobil))


x<-seq(-10, 42, by = 1)
length(x)
df3$observation<-NA
df3$observation[df3$type=="Above Mean"] <- x
df3$observation[df3$type=="Below Mean"] <- x

df3_covid$observation<-NA
df3_covid$observation[df3_covid$type=="Above Mean"] <- x
df3_covid$observation[df3_covid$type=="Below Mean"] <- x



df3$new_type<-NA
df3$new_type[as.numeric(df3$week)<11 & df3$type == "Above Mean"]<-"Above Mean Pre"
df3$new_type[as.numeric(df3$week)<11 & df3$type == "Below Mean"]<-"Below Mean Pre"
df3$new_type[as.numeric(df3$week)>=11 & df3$type == "Above Mean"]<-"Above Mean Post"
df3$new_type[as.numeric(df3$week)>=11 & df3$type == "Below Mean"]<-"Below Mean Post"


df3_covid$new_type<-NA
df3_covid$new_type[as.numeric(df3_covid$week)<11 & df3_covid$type == "Above Mean"]<-"Above Mean Pre"
df3_covid$new_type[as.numeric(df3_covid$week)<11 & df3_covid$type == "Below Mean"]<-"Below Mean Pre"
df3_covid$new_type[as.numeric(df3_covid$week)>=11 & df3_covid$type == "Above Mean"]<-"Above Mean Post"
df3_covid$new_type[as.numeric(df3_covid$week)>=11 & df3_covid$type == "Below Mean"]<-"Below Mean Post"



df3$week
as.numeric(df3$week)
df3_small<-subset(df3, as.numeric(df3$week)>=5 & as.numeric(df3$week)<=16)
irislabs1<-seq(min(as.numeric(df3_small$week)), max(as.numeric(df3_small$week),1))


df3$week<-as.numeric(df3$week)
sample_abom_bel10<-subset(df3, week<=10 & type=="Above Mean")
sample_abom_abo10<-subset(df3, week>10 & type=="Above Mean")
sample_belm_bel10<-subset(df3, week<=10 & type=="Below Mean")
sample_belm_abo10<-subset(df3, week>10 & type=="Below Mean")


df3_covid$week<-as.numeric(df3_covid$week)
sample_abom_bel10_covid<-subset(df3_covid, week<=10 & type=="Above Mean")
sample_abom_abo10_covid<-subset(df3_covid, week>=10 & type=="Above Mean")
sample_belm_bel10_covid<-subset(df3_covid, week<=10 & type=="Below Mean")
sample_belm_abo10_covid<-subset(df3_covid, week>=10 & type=="Below Mean")

mo_abom_bel10 <- lm(mean_mobil ~ week, data=sample_abom_bel10)
summary(mo_abom_bel10)
xseq<-seq(from=min(sample_abom_bel10$week), to=max(sample_abom_bel10$week))
pred <- predict(mo_abom_bel10)
df_abom_bel10 <- data.frame(week = xseq, pred)
df_abom_bel10$type_specific<-"Above Mean Pre"
df_abom_bel10$type<-"Above Mean"

mo_abom_abo10 <- lm(mean_mobil ~ week, data=sample_abom_abo10)
xseq<-seq(from=min(sample_abom_abo10$week), to=max(sample_abom_abo10$week))
pred <- predict(mo_abom_abo10)
df_abom_abo10 <- data.frame(week = xseq, pred)
df_abom_abo10$type_specific<-"Above Mean Post"
df_abom_abo10$type<-"Above Mean"


mo_belm_bel10 <- lm(mean_mobil ~ week, data=sample_belm_bel10)
xseq<-seq(from=min(sample_belm_bel10$week), to=max(sample_belm_bel10$week))
pred <- predict(mo_belm_bel10)
df_belm_bel10 <- data.frame(week = xseq, pred)
df_belm_bel10$type_specific<-"Below Mean Pre"
df_belm_bel10$type<-"Below Mean"

mo_belm_abo10 <- lm(mean_mobil ~ week, data=sample_belm_abo10)
xseq<-seq(from=min(sample_belm_abo10$week), to=max(sample_belm_abo10$week))
pred <- predict(mo_belm_abo10)
df_belm_abo10 <- data.frame(week = xseq, pred)
df_belm_abo10$type_specific<-"Below Mean Post"
df_belm_abo10$type<-"Below Mean"


new<-rbind(df_abom_bel10,
           df_abom_abo10,
           df_belm_bel10,
           df_belm_abo10)
new$line<-"predicted"
df3$line<-"actual"
names(new)[names(new)=="pred"]<-"mean_mobil"
names(new)[names(new)=="type_specific"]<-"new_type"


df3 <- df3[, c("week", "mean_mobil", "new_type", "type", "line")]
head(df3, 3)
head(new, 3)
df4<-rbind(df3, new)
x<-seq(-10, 42, by = 1)
length(x)
df4$observation<-NA
df4$observation[df4$type=="Above Mean"] <- x
df4$observation[df4$type=="Below Mean"] <- x


df4$observation<-as.character(df4$observation)
df4_covid$observation<-as.character(df4_covid$observation)

selection<-df4$observation[!stringr::str_detect(df4$observation, '-')]
selection_covid<-df4_covid$observation[!stringr::str_detect(df4_covid$observation, '-')]


df4$observation[!stringr::str_detect(df4$observation, '-')]<-paste("+", selection, sep = "")
df4$t<-paste("t", df4$observation, sep = "")
df4$t[df4$t=="t+0"]<-"t"


irislabs1<-seq(1,53,6)
irislabs2<-unique(df4[as.numeric(df4$week) %in% irislabs1,]$t)
length(irislabs2)
length(irislabs1)

df4_a<-df4
df4_a$mean_mobil[df4_a$type=="Below Mean"]<-NA
df4_a$mean_mobil[df4_a$week>11]<-NA
df4_a$mean_mobil[df4_a$line=="predicted"]<-NA


####################
#Trend All bonding#
####################

trend_all_bon<-ggplot(df4, aes(x=as.numeric(week), y=mean_mobil, linetype = line, color=type))+
  geom_line()+
  scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
trend_all_bon<-reposition_legend(trend_all_bon, 'bottom left')
ggsave(trend_all_bon, file = "./paper/graphs/figure_a7c.jpg", 
       height = 12, width = 20, 
       units = "cm", dpi = 200)

####################
#Trend Zoom bonding#
####################
df3_small<-subset(df4, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df3_small[as.numeric(df3_small$week) %in% irislabs1,]$t)
length(irislabs2)


trend_zoom_bon<-ggplot(df3_small, aes(x=as.numeric(week), y=mean_mobil, linetype = line, color=type))+
  geom_line()+
  scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
trend_zoom_bon<-reposition_legend(trend_zoom_bon, 'bottom left')
ggsave(trend_zoom_bon, file = "./paper/graphs/figure_a7a.jpg", 
       height = 12, width = 20, 
       units = "cm", dpi = 200)


#####################
#Vereine Time Trends#
#####################

df2 <- subset(data, select = c(mean_mobil, week, 
                               vereine_tot_pop_dummy))
df2<-subset(df2, !is.na(vereine_tot_pop_dummy))
df2$type<-ifelse(df2$vereine_tot_pop_dummy==1, "Above Mean", "Below Mean")
unique(df2$type)

df2<-st_drop_geometry(df2)
df3 <- df2%>%
  group_by(type, week)%>%
  summarise(mean_mobil = mean(mean_mobil))


x<-seq(-10, 42, by = 1)
length(x)
df3$observation<-NA
df3$observation[df3$type=="Above Mean"] <- x
df3$observation[df3$type=="Below Mean"] <- x

df3_covid$observation<-NA
df3_covid$observation[df3_covid$type=="Above Mean"] <- x
df3_covid$observation[df3_covid$type=="Below Mean"] <- x



df3$new_type<-NA
df3$new_type[as.numeric(df3$week)<11 & df3$type == "Above Mean"]<-"Above Mean Pre"
df3$new_type[as.numeric(df3$week)<11 & df3$type == "Below Mean"]<-"Below Mean Pre"
df3$new_type[as.numeric(df3$week)>=11 & df3$type == "Above Mean"]<-"Above Mean Post"
df3$new_type[as.numeric(df3$week)>=11 & df3$type == "Below Mean"]<-"Below Mean Post"


df3_covid$new_type<-NA
df3_covid$new_type[as.numeric(df3_covid$week)<11 & df3_covid$type == "Above Mean"]<-"Above Mean Pre"
df3_covid$new_type[as.numeric(df3_covid$week)<11 & df3_covid$type == "Below Mean"]<-"Below Mean Pre"
df3_covid$new_type[as.numeric(df3_covid$week)>=11 & df3_covid$type == "Above Mean"]<-"Above Mean Post"
df3_covid$new_type[as.numeric(df3_covid$week)>=11 & df3_covid$type == "Below Mean"]<-"Below Mean Post"



df3$week
as.numeric(df3$week)
df3_small<-subset(df3, as.numeric(df3$week)>=5 & as.numeric(df3$week)<=16)
irislabs1<-seq(min(as.numeric(df3_small$week)), max(as.numeric(df3_small$week),1))


df3$week<-as.numeric(df3$week)
sample_abom_bel10<-subset(df3, week<=10 & type=="Above Mean")
sample_abom_abo10<-subset(df3, week>10 & type=="Above Mean")
sample_belm_bel10<-subset(df3, week<=10 & type=="Below Mean")
sample_belm_abo10<-subset(df3, week>10 & type=="Below Mean")


df3_covid$week<-as.numeric(df3_covid$week)
sample_abom_bel10_covid<-subset(df3_covid, week<=10 & type=="Above Mean")
sample_abom_abo10_covid<-subset(df3_covid, week>=10 & type=="Above Mean")
sample_belm_bel10_covid<-subset(df3_covid, week<=10 & type=="Below Mean")
sample_belm_abo10_covid<-subset(df3_covid, week>=10 & type=="Below Mean")

mo_abom_bel10 <- lm(mean_mobil ~ week, data=sample_abom_bel10)
summary(mo_abom_bel10)
xseq<-seq(from=min(sample_abom_bel10$week), to=max(sample_abom_bel10$week))
pred <- predict(mo_abom_bel10)
df_abom_bel10 <- data.frame(week = xseq, pred)
df_abom_bel10$type_specific<-"Above Mean Pre"
df_abom_bel10$type<-"Above Mean"

mo_abom_abo10 <- lm(mean_mobil ~ week, data=sample_abom_abo10)
xseq<-seq(from=min(sample_abom_abo10$week), to=max(sample_abom_abo10$week))
pred <- predict(mo_abom_abo10)
df_abom_abo10 <- data.frame(week = xseq, pred)
df_abom_abo10$type_specific<-"Above Mean Post"
df_abom_abo10$type<-"Above Mean"


mo_belm_bel10 <- lm(mean_mobil ~ week, data=sample_belm_bel10)
xseq<-seq(from=min(sample_belm_bel10$week), to=max(sample_belm_bel10$week))
pred <- predict(mo_belm_bel10)
df_belm_bel10 <- data.frame(week = xseq, pred)
df_belm_bel10$type_specific<-"Below Mean Pre"
df_belm_bel10$type<-"Below Mean"

mo_belm_abo10 <- lm(mean_mobil ~ week, data=sample_belm_abo10)
xseq<-seq(from=min(sample_belm_abo10$week), to=max(sample_belm_abo10$week))
pred <- predict(mo_belm_abo10)
df_belm_abo10 <- data.frame(week = xseq, pred)
df_belm_abo10$type_specific<-"Below Mean Post"
df_belm_abo10$type<-"Below Mean"


new<-rbind(df_abom_bel10,
           df_abom_abo10,
           df_belm_bel10,
           df_belm_abo10)
new$line<-"predicted"
df3$line<-"actual"
names(new)[names(new)=="pred"]<-"mean_mobil"
names(new)[names(new)=="type_specific"]<-"new_type"


df3 <- df3[, c("week", "mean_mobil", "new_type", "type", "line")]
head(df3, 3)
head(new, 3)
df4<-rbind(df3, new)
x<-seq(-10, 42, by = 1)
length(x)
df4$observation<-NA
df4$observation[df4$type=="Above Mean"] <- x
df4$observation[df4$type=="Below Mean"] <- x


df4$observation<-as.character(df4$observation)
df4_covid$observation<-as.character(df4_covid$observation)

selection<-df4$observation[!stringr::str_detect(df4$observation, '-')]
selection_covid<-df4_covid$observation[!stringr::str_detect(df4_covid$observation, '-')]


df4$observation[!stringr::str_detect(df4$observation, '-')]<-paste("+", selection, sep = "")
df4$t<-paste("t", df4$observation, sep = "")
df4$t[df4$t=="t+0"]<-"t"


irislabs1<-seq(1,53,6)
irislabs2<-unique(df4[as.numeric(df4$week) %in% irislabs1,]$t)
length(irislabs2)
length(irislabs1)

df4_a<-df4
df4_a$mean_mobil[df4_a$type=="Below Mean"]<-NA
df4_a$mean_mobil[df4_a$week>11]<-NA
df4_a$mean_mobil[df4_a$line=="predicted"]<-NA


####################
#Trend All vereine#
####################

trend_all_ver<-ggplot(df4, aes(x=as.numeric(week), y=mean_mobil, linetype = line, color=type))+
  geom_line()+
  scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
trend_all_ver<-reposition_legend(trend_all_ver, 'bottom left')
ggsave(trend_all_ver, file = "./paper/graphs/figure_a8c.jpg", 
       height = 12, width = 20, 
       units = "cm", dpi = 200)

####################
#Trend Zoom vereine#
####################
df3_small<-subset(df4, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df3_small[as.numeric(df3_small$week) %in% irislabs1,]$t)
length(irislabs2)

trend_zoom_ver<-ggplot(df3_small, aes(x=as.numeric(week), y=mean_mobil, linetype = line, color=type))+
  geom_line()+
  scale_linetype_manual(name= "", values = c("predicted"="dashed", "actual"="solid")) +
  guides(linetype = "none")+
  scale_colour_manual(name="", values=c("Above Mean"="red", "Below Mean"="blue"))+
  geom_vline(xintercept = 10, color ="red")+
  #scale_x_continuous(name = "Weeks", breaks=seq(1,55,5))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  labs(x = "Week", y = "Mobility Compared to 2019")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
trend_zoom_ver<-reposition_legend(trend_zoom_ver, 'bottom left')
ggsave(trend_zoom_ver, file = "./paper/graphs/figure_a8a.jpg", 
       height = 12, width = 20, 
       units = "cm", dpi = 200)



###########
#All clubs#
###########

#Create dummies
library('fastDummies')
# Create dummy variables:
dataf <- dummy_cols(data, select_columns = 'AGS')
dataw <- dummy_cols(dataf, select_columns = 'week')
datas <- dummy_cols(dataw, select_columns = 'SN_L')

library(dplyr)
new<-dataw %>% dplyr::select(starts_with("AGS_"))
new2<-dataw %>% dplyr::select(starts_with("week_"))
new3<-datas %>% dplyr::select(starts_with("SN_L"))

list_dummies<-names(new)
#I am removing the counties which get dropped in Stata

#Re-enable this later if necessary
list_dummies<-list_dummies[list_dummies!= "AGS_0"]
list_dummies<-list_dummies[list_dummies!= "AGS_5316"]
list_dummies<-list_dummies[list_dummies!= "AGS_8121"]
list_dummies<-list_dummies[list_dummies!= "AGS_8221"]
list_dummies<-list_dummies[list_dummies!= "AGS_8231"]
list_dummies<-list_dummies[list_dummies!= "AGS_8311"]
list_dummies<-list_dummies[list_dummies!= "AGS_8421"]
list_dummies<-list_dummies[list_dummies!= "AGS_9461"]
list_dummies<-list_dummies[list_dummies!= "AGS_13004"]
list_dummies<-list_dummies[list_dummies!= "AGS_16052"]
list_dummies<-list_dummies[list_dummies!= "AGS_16053"]
list_dummies<-list_dummies[list_dummies!= "AGS_16054"]
list_dummies<-list_dummies[list_dummies!= "AGS_16055"]
list_dummies<-list_dummies[list_dummies!= "AGS_16061"]
list_dummies<-list_dummies[list_dummies!= "AGS_16077"]

list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]

list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_09"]
list_week_dummies2<-list_week_dummies2[list_week_dummies2!= "week_52"]


list_land_dummies<-names(new3)
list_land_dummies_full<-list_land_dummies[list_land_dummies!= "SN_L"]
list_land_dummies<-list_land_dummies_full[list_land_dummies_full!= "SN_L_01"]

#Creating the interaction terms for bonding
for (i in list_land_dummies){
  for (j in list_week_dummies2)
    datas[[paste("inter_", i, "_", j, sep="")]]<-datas[[i]]*datas[[j]]
}

new<-datas %>% dplyr::select(starts_with("inter_"))
list_land_dummies_star<-names(new)

#Event study
#Bridging Interaction
new2<-dataw %>% dplyr::select(starts_with("week_"))
list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]

#Creating the interaction terms for bridging
for (i in list_week_dummies){
  datas[[paste("bri_", i, sep="")]]<-datas[[i]]*datas$bridging_tot_pop_dummy
}

#Creating the interaction terms for bonding
for (i in list_week_dummies){
  datas[[paste("bon_", i, sep="")]]<-datas[[i]]*datas$bonding_tot_pop_dummy
}

#Creating the interaction terms for all vereine
for (i in list_week_dummies){
  datas[[paste("ver_", i, sep="")]]<-datas[[i]]*datas$vereine_tot_pop_dummy
}


#Bri
inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_week_10"]

#Bon
inter_terms_bon<-datas %>% dplyr::select(starts_with("bon_week"))
inter_terms_bon<-names(inter_terms_bon)
inter_terms_bon<-inter_terms_bon[inter_terms_bon != "bon_week_10"]


#Ver
inter_terms_ver<-datas %>% dplyr::select(starts_with("ver_week"))
inter_terms_ver<-names(inter_terms_ver)
inter_terms_ver<-inter_terms_ver[inter_terms_ver != "ver_week_10"]


#Deleting variables which are collinear
list_land_dummies_full2<-list_land_dummies_full[list_land_dummies_full!= "SN_L_02"]
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_53"]

for_event_bri<-list.append(inter_terms_bri, list_land_dummies_full, list_week_dummies)
for_event_bri2<-list.append(inter_terms_bri, list_land_dummies_full2, list_week_dummies2)

for_event_bri_covid<-list.append("covid_diff", inter_terms_bri, list_land_dummies_full, list_week_dummies)
for_event_bri2_covid<-list.append("covid_diff", inter_terms_bri, list_land_dummies_full2, list_week_dummies2)

for_event_bon<-list.append(inter_terms_bon, list_land_dummies_full, list_week_dummies)
for_event_bon2<-list.append(inter_terms_bon, list_land_dummies_full2, list_week_dummies2)

for_event_ver<-list.append(inter_terms_ver, list_land_dummies_full, list_week_dummies)
for_event_ver2<-list.append(inter_terms_ver, list_land_dummies_full2, list_week_dummies2)

#Deleting missing observations
datas$time<-as.numeric(datas$week)
case_bri<-subset(datas, select = list.append(for_event_bri2, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y",
                                         "time", "AGS"))
case_bri<-case_bri[complete.cases(case_bri), ]

case_bri_covid<-subset(datas, select = list.append(for_event_bri2_covid, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y",
                                             "time", "AGS"))
case_bri_covid<-case_bri_covid[complete.cases(case_bri_covid), ]

summary(lm(covid_diff~bridging_total, data=case_bri_covid))
summary(lm(mean_mobil~bridging_total, data=case_bri))


case_bon<-subset(datas, select = list.append(for_event_bon2, "bonding_total", "mean_mobil", "POINT_X", "POINT_Y",
                                             "time", "AGS"))
case_bon<-case_bon[complete.cases(case_bon), ]
case_ver<-subset(datas, select = list.append(for_event_ver2, "gvereine", "mean_mobil", "POINT_X", "POINT_Y",
                                             "time", "AGS"))
case_ver<-case_ver[complete.cases(case_ver), ]

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(case_bri$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Merging with spatial dataframe
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(as.character(nuts2$AGS))
case_bri_oneweek<-subset(case_bri, week_11==1)
case_bri_oneweek_covid<-subset(case_bri_covid, week_11==1)

case_bon_oneweek<-subset(case_bon, week_11==1)
case_ver_oneweek<-subset(case_ver, week_11==1)

#Bridging
nuts_merged<-left_join(nuts2, case_bri_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

#Bridging Covid
nuts_merged<-left_join(nuts2, case_bri_oneweek_covid, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]


#Bonding
nuts_merged<-left_join(nuts2, case_bon_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bonding_total))

#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]


#Vereine
nuts_merged<-left_join(nuts2, case_ver_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$gvereine.x))
nuts_merged2$gvereine<-as.numeric(nuts_merged2$gvereine.x)

#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

#Getting coordinates
coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)


#Bri formula
for_event_bri<-list.append(inter_terms_bri, list_land_dummies_full, list_week_dummies)
time_trend_bri<-as.formula(paste("mean_mobil~",
                                 paste(for_event_bri2, collapse= "+")))


time_trend_bri_covid<-as.formula(paste("covid_diff~",
                                 paste(for_event_bri2, collapse= "+")))



#Bon formula
for_event_bon<-list.append(inter_terms_bon, list_land_dummies_full, list_week_dummies)
time_trend_bon<-as.formula(paste("mean_mobil~",
                                 paste(for_event_bon2, collapse= "+")))

#Ver formula
for_event_ver<-list.append(inter_terms_ver, list_land_dummies_full, list_week_dummies)
time_trend_ver<-as.formula(paste("mean_mobil~",
                                 paste(for_event_ver2, collapse= "+")))


#BRIDGING: Creating distance matrix - accelerates the conely regression
dm_bri <- dist_mat(case_bri, lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time")

#Takes 30 sec
df_bri<-conleyreg(time_trend_bri, data=case_bri, 455, unit = "AGS", time = "time", 
                  dist_mat=dm_bri, lag_cutoff=temp_lag)
df_bri<-broom::tidy(df_bri)
df_bri$term[grepl( "bri_week_" , df_bri$term)]
df_bri<-add_row(df_bri,term="bri_week_10",estimate=0)
df_bri<-subset(df_bri, grepl( "bri_week_" , df_bri$term))
df_bri<-df_bri[order(df_bri$term),]

x<-seq(-10, 42, by = 1)
length(x)
df_bri$observation<-NA
df_bri$observation <- x

df_bri$observation<-as.character(df_bri$observation)
selection<-df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]
df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]<-paste("+", selection, sep = "")
df_bri$t<-paste("Network x \n", "(t", df_bri$observation, ")", sep = "")
df_bri$t[df_bri$t=="Network x \n(t+0)"]<-"Network x t"

df_bri$week<-NA
df_bri$week<-readr::parse_number(df_bri$term)
df_bri$model<-"Bridging"


#Takes 30 sec
summary(case_bri_covid$covid_diff)
df_bri_covid<-conleyreg(time_trend_bri_covid, data=case_bri_covid, 455, unit = "AGS", time = "time", 
                  dist_mat=dm_bri, lag_cutoff=temp_lag)
df_bri_covid<-broom::tidy(df_bri_covid)
df_bri_covid$term[grepl( "bri_week_" , df_bri_covid$term)]
df_bri_covid<-add_row(df_bri_covid,term="bri_week_10",estimate=0)
df_bri_covid<-subset(df_bri_covid, grepl( "bri_week_" , df_bri_covid$term))
df_bri_covid<-df_bri_covid[order(df_bri_covid$term),]

x<-seq(-10, 42, by = 1)
length(x)
df_bri_covid$observation<-NA
df_bri_covid$observation <- x

df_bri_covid$observation<-as.character(df_bri_covid$observation)
selection<-df_bri_covid$observation[!stringr::str_detect(df_bri_covid$observation, '-')]
df_bri_covid$observation[!stringr::str_detect(df_bri_covid$observation, '-')]<-paste("+", selection, sep = "")
df_bri_covid$t<-paste("Network x \n", "(t", df_bri_covid$observation, ")", sep = "")
df_bri_covid$t[df_bri_covid$t=="Network x \n(t+0)"]<-"Network x t"

df_bri_covid$week<-NA
df_bri_covid$week<-readr::parse_number(df_bri_covid$term)
df_bri_covid$model<-"Bridging"


#Bonding: Creating distance matrix - accelerates the conely regression
dm_bon <- dist_mat(case_bon, lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time")
#Takes 30 sec
df_bon<-conleyreg(time_trend_bon, data=case_bon, 538, unit = "AGS", time = "time", 
                  dist_mat=dm_bon, lag_cutoff=temp_lag)
df_bon<-broom::tidy(df_bon)
df_bon$term[grepl( "bon_week_" , df_bon$term)]
df_bon<-add_row(df_bon,term="bon_week_10",estimate=0)
df_bon<-subset(df_bon, grepl( "bon_week_" , df_bon$term))
df_bon<-df_bon[order(df_bon$term),]

x<-seq(-10, 42, by = 1)
length(x)
df_bon$observation<-NA
df_bon$observation <- x

df_bon$observation<-as.character(df_bon$observation)
selection<-df_bon$observation[!stringr::str_detect(df_bon$observation, '-')]
df_bon$observation[!stringr::str_detect(df_bon$observation, '-')]<-paste("+", selection, sep = "")
df_bon$t<-paste("Network x \n", "(t", df_bon$observation, ")", sep = "")
df_bon$t[df_bon$t=="Network x \n(t+0)"]<-"Network x t"

df_bon$week<-NA
df_bon$week<-readr::parse_number(df_bon$term)
df_bon$model<-"Bonding"

#Vereine: Creating distance matrix - accelerates the conely regression
dm_ver <- dist_mat(case_ver, lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time")
#Takes 30 sec
df_ver<-conleyreg(time_trend_ver, data=case_ver, 455, unit = "AGS", time = "time", 
                  dist_mat=dm_ver, lag_cutoff=temp_lag)
df_ver<-broom::tidy(df_ver)
df_ver$term[grepl( "ver_week_" , df_ver$term)]
df_ver<-add_row(df_ver,term="ver_week_10",estimate=0)
df_ver<-subset(df_ver, grepl( "ver_week_" , df_ver$term))
df_ver<-df_ver[order(df_ver$term),]

x<-seq(-10, 42, by = 1)
length(x)
df_ver$observation<-NA
df_ver$observation <- x

df_ver$observation<-as.character(df_ver$observation)
selection<-df_ver$observation[!stringr::str_detect(df_ver$observation, '-')]
df_ver$observation[!stringr::str_detect(df_ver$observation, '-')]<-paste("+", selection, sep = "")
df_ver$t<-paste("Network x \n", "(t", df_ver$observation, ")", sep = "")
df_ver$t[df_ver$t=="Network x \n(t+0)"]<-"Network x t"

df_ver$week<-NA
df_ver$week<-readr::parse_number(df_ver$term)
df_ver$model<-"Vereine"

##########################
#All Event Study Bridging#
##########################
irislabs1<-seq(1,53,6)
irislabs2<-unique(df_bri[as.numeric(df_bri$week) %in% irislabs1,]$t)
length(irislabs2)

coef_all_bri<-ggplot(df_bri, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(coef_all_bri, file = "./paper/graphs/figure_1d.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


##########################
#All Event Study Bridging#
##########################
irislabs1<-seq(5,16,1)
irislabs2<-unique(df_bri[as.numeric(df_bri$week) %in% irislabs1,]$t)
length(irislabs2)
df_bri_zoom<-subset(df_bri, week>4 & week<17)

coef_zoom_bri<-ggplot(df_bri_zoom, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(coef_zoom_bri, file = "./paper/graphs/figure_1b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



################################
#All Event Study Bridging COVID#
################################
irislabs1<-seq(1,53,6)
irislabs2<-unique(df_bri[as.numeric(df_bri_covid$week) %in% irislabs1,]$t)
length(irislabs2)

coef_all_bri_covid<-ggplot(df_bri_covid, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(coef_all_bri_covid, file = "./paper/graphs/figure_a6b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



###########################
#Zoom Event Study Bridging#
###########################

df_bri_small<-subset(df_bri, week>=(5) & week<=(16))

irislabs1<-seq(1,53,1)
irislabs2<-unique(df_bri[as.numeric(df_bri$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_bri<-ggplot(df_bri_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(coef_zoom_bri, file = "./paper/graphs/figure_a9a.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


##########################
#All Event Study Bonding#
##########################
irislabs1<-seq(1,53,6)
irislabs2<-unique(df_bon[as.numeric(df_bon$week) %in% irislabs1,]$t)
length(irislabs2)

coef_all_bon<-ggplot(df_bon, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(coef_all_bon, file = "./paper/graphs/figure_a7d.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



###########################
#Zoom Event Study Bonding#
###########################

df_bon_small<-subset(df_bon, week>=(5) & week<=(16))

irislabs1<-seq(1,53,1)
irislabs2<-unique(df_bon[as.numeric(df_bon$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_bon<-ggplot(df_bon_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_bon, file = "./paper/graphs/figure_a7b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



##########################
#All Event Study Vereine#
##########################
irislabs1<-seq(1,53,6)
irislabs2<-unique(df_ver[as.numeric(df_ver$week) %in% irislabs1,]$t)
length(irislabs2)

coef_all_ver<-ggplot(df_ver, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(coef_all_ver, file = "./paper/graphs/figure_a8d.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)



###########################
#Zoom Event Study Vereine#
###########################

df_ver_small<-subset(df_ver, week>=(5) & week<=(16))

irislabs1<-seq(1,53,1)
irislabs2<-unique(df_ver[as.numeric(df_ver$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_ver<-ggplot(df_ver_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
ggsave(coef_zoom_ver, file = "./paper/graphs/figure_a8b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)

